"""
Mathematical validation tests for vector and tensor fields

Tests field operations against analytical solutions to ensure correctness.
These tests verify that the numerical implementations produce the correct
results for known mathematical identities and analytical solutions.
"""

using Test
using LinearAlgebra
using PDEJulia
using PDEJulia.Grids
using PDEJulia.Fields

@testset "Mathematical Validation Tests" begin

@testset "Vector Field Analytical Validation" begin
    # Use a finer grid for better accuracy in derivative tests
    grid_2d = UnitGrid([41, 41])
    h = spacing(grid_2d)[1]  # Grid spacing
    
    @testset "Divergence Validation" begin
        # Test div(grad(φ)) = ∇²φ for scalar field φ = x² + y²
        # grad(φ) = [2x, 2y]
        # div(grad(φ)) = 2 + 2 = 4 everywhere
        
        scalar_phi = ScalarField(pos -> pos[1]^2 + pos[2]^2, grid_2d)
        grad_phi = gradient(scalar_phi)
        div_grad_phi = divergence(grad_phi)
        
        # Should be 4 everywhere (except near boundaries due to finite differences)
        interior_mask = 3:(size(grid_2d)[1]-2), 3:(size(grid_2d)[2]-2)  # Avoid boundary effects
        interior_values = data(div_grad_phi)[interior_mask...]
        
        @test all(abs.(interior_values .- 4.0) .< 0.1)  # Allow small numerical error
    end
    
    @testset "Curl Validation" begin
        # Test curl of conservative field should be zero
        # For φ = xy, grad(φ) = [y, x] is conservative
        # curl([y, x]) = ∂x/∂x - ∂y/∂y = 1 - 1 = 0
        
        conservative_field = VectorField(pos -> [pos[2], pos[1]], grid_2d)
        curl_result = curl(conservative_field)
        
        # Should be zero everywhere (except boundary effects)
        interior_mask = 3:(size(grid_2d)[1]-2), 3:(size(grid_2d)[2]-2)
        interior_values = data(curl_result)[interior_mask...]
        
        @test all(abs.(interior_values) .< 0.1)
    end
    
    @testset "Gradient-Divergence Identity" begin
        # Test vector identity: div(f*v) = f*div(v) + grad(f)·v
        # where f is scalar field and v is vector field
        
        f = ScalarField(pos -> sin(pos[1]) * cos(pos[2]), grid_2d)
        v = VectorField(pos -> [pos[1], pos[2]^2], grid_2d)
        
        # Left side: div(f*v)
        fv = f * v
        div_fv = divergence(fv)
        
        # Right side: f*div(v) + grad(f)·v
        div_v = divergence(v)
        f_div_v = f * div_v
        
        grad_f = gradient(f)
        grad_f_dot_v = dot(grad_f, v)
        
        right_side = f_div_v + grad_f_dot_v
        
        # Compare interior points
        interior_mask = 3:(size(grid_2d)[1]-2), 3:(size(grid_2d)[2]-2)
        left_interior = data(div_fv)[interior_mask...]
        right_interior = data(right_side)[interior_mask...]
        
        @test all(abs.(left_interior - right_interior) .< 0.2)
    end
    
    @testset "Vector Cross Product Properties" begin
        grid_3d = UnitGrid([21, 21, 21])
        
        # Test anti-commutativity: a × b = -(b × a)
        a = VectorField(pos -> [1.0, pos[1], pos[2]], grid_3d)
        b = VectorField(pos -> [pos[3], 1.0, pos[1]], grid_3d)
        
        cross_ab = cross(a, b)
        cross_ba = cross(b, a)
        neg_cross_ba = -cross_ba
        
        @test all(abs.(data(cross_ab) - data(neg_cross_ba)) .< 1e-14)
        
        # Test that a × a = 0
        cross_aa = cross(a, a)
        @test all(abs.(data(cross_aa)) .< 1e-14)
    end
end

@testset "Tensor Field Analytical Validation" begin
    grid_2d = UnitGrid([21, 21])
    
    @testset "Tensor Invariants Validation" begin
        # Create a known 2x2 matrix field: T = [a b; c d]
        # with known eigenvalues and invariants
        
        tensor_field = Tensor2Field(pos -> [2.0 1.0; 0.0 3.0], grid_2d)  # Upper triangular
        
        # Test trace (first invariant)
        trace_result = trace(tensor_field)
        @test all(data(trace_result) .≈ 5.0)  # 2 + 3 = 5
        
        # Test determinant (second invariant in 2D)
        det_result = det(tensor_field)
        @test all(data(det_result) .≈ 6.0)  # 2*3 - 1*0 = 6
        
        # For upper triangular matrix, eigenvalues are diagonal elements
        eigenvals_fields = eigenvals(tensor_field)
        @test all(data(eigenvals_fields[1]) .≈ 3.0)  # Larger eigenvalue
        @test all(data(eigenvals_fields[2]) .≈ 2.0)  # Smaller eigenvalue
    end
    
    @testset "Tensor Decomposition" begin
        # Test symmetric-antisymmetric decomposition
        # Any tensor T can be written as T = S + A where S = (T + T^T)/2, A = (T - T^T)/2
        
        tensor_field = Tensor2Field(pos -> [1.0 3.0; 2.0 4.0], grid_2d)
        
        # Symmetric part
        symmetric_part = symmetrize(tensor_field)
        
        # Antisymmetric part  
        antisymmetric_part = tensor_field - symmetric_part
        
        # Check that symmetric_part is actually symmetric
        @test is_symmetric(symmetric_part)
        
        # Check that antisymmetric_part has zero trace
        antisym_trace = trace(antisymmetric_part)
        @test all(abs.(data(antisym_trace)) .< 1e-14)
        
        # Check that sum recovers original
        reconstructed = symmetric_part + antisymmetric_part
        @test all(abs.(data(reconstructed) - data(tensor_field)) .< 1e-14)
    end
    
    @testset "Tensor-Vector Product Validation" begin
        # Test (AB)v = A(Bv) for tensors A, B and vector v
        grid_2d_small = UnitGrid([11, 11])  # Smaller grid for efficiency
        
        A = Tensor2Field(pos -> [2.0 1.0; 0.0 3.0], grid_2d_small)
        B = Tensor2Field(pos -> [1.0 2.0; 1.0 1.0], grid_2d_small)
        v = VectorField(pos -> [1.0, 2.0], grid_2d_small)
        
        # Left side: (AB)v
        # Note: We need tensor-tensor multiplication, which we'll implement as component-wise
        AB_v = dot(A, dot(B, v))
        
        # Right side: A(Bv)  
        Bv = dot(B, v)
        A_Bv = dot(A, Bv)
        
        @test all(abs.(data(AB_v) - data(A_Bv)) .< 1e-12)
    end
end

@testset "Field Collection Validation" begin
    grid_2d = UnitGrid([21, 21])
    
    @testset "Conservation Properties" begin
        # Test that field collection operations preserve total mass/energy
        
        # Create reaction-diffusion system
        fields = reaction_diffusion_fields(grid_2d, 2)
        
        # Initialize with some mass
        fields[1] = ScalarField(pos -> exp(-(pos[1]^2 + pos[2]^2)), grid_2d)
        fields[2] = ScalarField(pos -> exp(-((pos[1]-0.5)^2 + (pos[2]-0.5)^2)), grid_2d)
        
        # Calculate initial total mass
        initial_integrals = integrate(fields)
        initial_total = sum(initial_integrals)
        
        # Apply some operations that should preserve mass
        scaled_fields = fields * 2.0
        scaled_integrals = integrate(scaled_fields)
        scaled_total = sum(scaled_integrals)
        
        @test abs(scaled_total - 2*initial_total) < 1e-12
        
        # Test that addition preserves structure
        added_fields = fields + fields
        added_integrals = integrate(added_fields)
        added_total = sum(added_integrals)
        
        @test abs(added_total - 2*initial_total) < 1e-12
    end
    
    @testset "Multi-Physics Coupling" begin
        # Test fluid-structure interaction-like coupling
        fluid_fields = fluid_dynamics_fields(grid_2d)
        elastic_fields = elasticity_fields(grid_2d)
        
        # Initialize with some test data
        fluid_fields["density"] = ScalarField(pos -> 1.0 + 0.1*sin(pos[1]), grid_2d)
        fluid_fields["velocity"] = VectorField(pos -> [pos[2], -pos[1]], grid_2d)
        
        elastic_fields["displacement"] = VectorField(pos -> [0.01*pos[1], 0.01*pos[2]], grid_2d)
        elastic_fields["stress"] = Tensor2Field(pos -> [1.0 0.0; 0.0 1.0], grid_2d)
        
        # Create combined system
        coupled_system = FieldCollection([
            fluid_fields["density"],
            fluid_fields["velocity"],
            elastic_fields["displacement"],
            elastic_fields["stress"]
        ], ["density", "velocity", "displacement", "stress"])
        
        @test length(coupled_system) == 4
        @test coupled_system["density"] isa ScalarField
        @test coupled_system["velocity"] isa VectorField
        @test coupled_system["displacement"] isa VectorField
        @test coupled_system["stress"] isa Tensor2Field
        
        # Test that operations work on mixed field types
        stats = get_statistics(coupled_system)
        @test stats.n_fields == 4
        @test length(unique(stats.field_types)) == 3  # Three different types
    end
end

@testset "Numerical Accuracy Validation" begin
    @testset "Convergence Testing" begin
        # Test that operations converge to analytical solutions as grid is refined
        
        grid_sizes = [11, 21, 41]
        errors = Float64[]
        
        for N in grid_sizes
            grid = UnitGrid([N, N])
            
            # Test gradient of quadratic function: ∇(x² + y²) = [2x, 2y]
            scalar_field = ScalarField(pos -> pos[1]^2 + pos[2]^2, grid)
            grad_field = gradient(scalar_field)
            
            # Analytical solution
            analytical_grad = VectorField(pos -> [2*pos[1], 2*pos[2]], grid)
            
            # Compute error (avoid boundaries where finite differences are less accurate)
            error = 0.0
            count = 0
            for i in 3:(N-2), j in 3:(N-2)
                numerical = [data(grad_field)[i, j, 1], data(grad_field)[i, j, 2]]
                analytical = [data(analytical_grad)[i, j, 1], data(analytical_grad)[i, j, 2]]
                error += norm(numerical - analytical)
                count += 1
            end
            error /= count
            push!(errors, error)
        end
        
        # Check that error decreases with finer grid (should be approximately quadratic)
        @test errors[2] < errors[1]  # 21x21 better than 11x11
        @test errors[3] < errors[2]  # 41x41 better than 21x21
        
        # Rough convergence rate check (finite differences should be ~O(h²))
        conv_rate_1 = log(errors[1]/errors[2]) / log(2.0)  # 11->21 is 2x refinement
        conv_rate_2 = log(errors[2]/errors[3]) / log(2.0)  # 21->41 is ~2x refinement
        
        # Should be approximately 2nd order
        @test 1.5 < conv_rate_1 < 3.0
        @test 1.5 < conv_rate_2 < 3.0
    end
    
    @testset "Stability Testing" begin
        # Test that operations remain stable and don't produce NaN or Inf
        
        grid_2d = UnitGrid([21, 21])
        
        # Test with extreme values
        large_scalar = ScalarField(pos -> 1e6, grid_2d)
        small_scalar = ScalarField(pos -> 1e-6, grid_2d)
        
        large_vector = VectorField(pos -> [1e6, 1e6], grid_2d)
        small_vector = VectorField(pos -> [1e-6, 1e-6], grid_2d)
        
        # Test operations don't produce NaN or Inf
        ops_to_test = [
            () -> large_scalar + small_scalar,
            () -> large_scalar * 1e-6,
            () -> dot(large_vector, small_vector),
            () -> magnitude(large_vector),
            () -> normalize(large_vector)
        ]
        
        for op in ops_to_test
            result = op()
            field_data = data(result)
            @test all(isfinite.(field_data))
            @test !any(isnan.(field_data))
        end
    end
end

@testset "Boundary Effects Testing" begin
    # Test that field operations handle boundaries appropriately
    
    grid_2d = UnitGrid([21, 21])
    
    @testset "Gradient Near Boundaries" begin
        # Create a field that has known gradient everywhere
        linear_field = ScalarField(pos -> 2*pos[1] + 3*pos[2], grid_2d)
        grad_field = gradient(linear_field)
        
        # Gradient should be [2, 3] everywhere, even near boundaries
        # (though numerical accuracy may be lower at boundaries)
        
        # Check interior points are accurate
        interior_x = data(grad_field)[5:17, 5:17, 1]
        interior_y = data(grad_field)[5:17, 5:17, 2]
        
        @test all(abs.(interior_x .- 2.0) .< 0.01)
        @test all(abs.(interior_y .- 3.0) .< 0.01)
        
        # Check that boundary points are reasonable (may have larger error)
        boundary_x = [data(grad_field)[1, :, 1]; data(grad_field)[end, :, 1]; 
                     data(grad_field)[:, 1, 1]; data(grad_field)[:, end, 1]]
        boundary_y = [data(grad_field)[1, :, 2]; data(grad_field)[end, :, 2]; 
                     data(grad_field)[:, 1, 2]; data(grad_field)[:, end, 2]]
        
        @test all(abs.(boundary_x .- 2.0) .< 1.0)  # Allow larger error at boundaries
        @test all(abs.(boundary_y .- 3.0) .< 1.0)
    end
end

end  # Main test set